
\chapter{数值积分}

我们前面说了，讲义第二部分是关于函数的一切。
这一章是关于函数的积分。
函数的积分方法可以分为两类；一类方法在积分区域内选取一系列等距的节点，在节点之间等长的小区间上做近似积分，然后加起来，这一类积分叫做“{\bf Newton-Cortes}"方法；
另一类方法在积分区域内选取一些特定的{\bf 不等距}的节点，然后用这些节点上的函数值，构造整个积分的近似值，这一类积分叫做{\bf 高斯积分}。

这两类积分各有优劣，对于震荡比较厉害的积分，应该用第一类，即Newton-Cortes方法；而对于比较平缓的积分，适用高斯积分，高斯积分可以在很少的资源下得到很高的精度。
在后文中，我们会写程序进行比较说明。

下面我们先说Newton-Cortes积分。

\section{Newton-Cortes积分}
如果需要求取积分
\begin{eqnarray}
I_f = \int^{b}_{a} f(x) dx,
\end{eqnarray}
其中 $I$ 是 $Integral$ 的缩写，就是求函数jizhi$f(x)$下方的面积，如图 \ref{figure:integration_demon}所示。
\begin{figure}[h]
\centering
    \includegraphics[width=0.8\textwidth]{integration_demon.eps}
    \caption{　f(x)积分示意图 }
    \label{figure:integration_demon}
\end{figure}

我们可以将积分区间 $\left[a, b\right]$ 分割为 $n$ 个均等的小区间
\begin{equation}
\left[x_0=a, x_1\right],\cdots,\left[x_{n-1}, x_n=b \right] \nonumber
\end{equation}
把整个积分分为很多小份
\begin{equation}
I_f = \sum^{n-1}_{i=0} \int^{x_{i+1}}_{x_i} f(x) dx.
\end{equation}

\paragraph{中值法}
我们可以将每个区间的积分近似为
\begin{equation}
\int^{x_{i+1}}_{x_i} f(x) dx \approx h f(\frac{x_i + x_{i+1}}{2}),
\end{equation}
当区间足够小的时候，小区间中的函数值都很接近，所以可以做上面的近似。
总积分则近似为
\begin{eqnarray}
I_f \approx \sum^{n-1}_{i=0} h f(\frac{x_i + x_{i+1}}{2}).
\end{eqnarray}
这个方法取每个小区间的中点的函数值，进行近似，所以叫做中值法。
\begin{lstlisting}
/*
 * integ_midpoint: 中值法积分
 * 给定区间[a,b]，分割为n份，被积函数为func，返回积分值
 */
double integ_midpoint(double a, double b, int n, double (*func)(double)){

        int i;
        double y=0;
        double h=(b-a)/n;
        for(i=0;i<n;i++){
                y += func( a + (i+0.5)*h );
        }
        y *= h;
        return y;
}
\end{lstlisting}
这里我省去了头文件、主函数、被积函数定义等内容，因为这些内容在前文中有示例。

\paragraph{梯形法}
我们也可以换一个方法做近似
\begin{equation}
\int^{x_{i+1}}_{x_i} f(x) dx \approx h (f(x_i) + f(x_{i+1}))/2,
\end{equation}
上式右边正是梯形面积公式。
这个近似的几何意义是：用一个梯形的面积对小区间内的积分做近似，所以这个方法叫做梯形法。
总的积分值为
\begin{eqnarray}
I_f &\approx& \sum^{n-1}_{i=0} h (f(x_i) + f(x_{i+1}))/2 \nonumber\\
    &=& h \{ \frac{f(x_0)}{2} + f(x_1) + \cdots + f(x_{n-1}) + \frac{f(x_n)}{2} \}.
\end{eqnarray}
\begin{lstlisting}
/*
 * integ_trapzoid: 梯形法积分
 * 给定区间[a,b]，分割为n份，被积函数为func，返回积分值
 */
double integ_trapzoid(double a, double b, int n, double (*func)(double)){

        int i;
        double y=func(a)/2;
        double h=(b-a)/n;
        for(i=1;i<=n-1;i++){
                y += func( a + i*h );
        }
        y += func(b)/2;
        y *= h;
        return y;
}
\end{lstlisting}

\paragraph{辛普生方法}
我们可以换一个角度看问题，
在小区间内，中值法是用一个常数作为函数的近似
\begin{equation}
f(x) \approx a,
\end{equation}
而梯形法是用１阶多项式作为函数的近似
\begin{equation}
f(x) \approx a + bx.
\end{equation}
那么，如果我们想要更精确的近似，可以取三个点
\begin{eqnarray}
&&(x_i, f(x_i)), \nonumber\\
&&(x_i+h/2, f(x_i+h/2)), \nonumber\\
&&(x_i+h, f(x_i+h))，
\end{eqnarray}
进行２阶多项式拟合。
利用前文中多项式拟合的拉格朗日公式，可以得到
\begin{eqnarray}
f(x) &\approx& f(x_i) \frac{ (x-x_{h/2})(x-x_h) }{ (-h/2)(-h) } \nonumber\\
    && + f(x_i + h/2) \frac{ (x-x_i)(x-x_h) }{ (h/2) (-h/2)} \nonumber\\
    && + f(x_i + h) \frac{ (x-x_i)(x-x_{h/2}) }{ h \cdot h/2},
\end{eqnarray}
为了方便，我们用$x_{h/2},x_h$代表$x_i + h/2, x_i + h$。
在区间 $\left[x_i, x_{i+1}\right]$ 内积分以后得到
\begin{equation}
\int^{x_i +h}_{x_i} f(x) dx \approx \frac{h}{6}[ f(x_i) + 4f(x_i + h/2) + f(x_i + h) ].
\end{equation}
所以有
\begin{eqnarray}
I_f &\approx& \sum^{n-1}_{i=0} \frac{h}{6} \left[ f(x_i) + 4 f(\frac{x_i + x_{i+1}}{2}) + f(x_{i+1}) \right] \nonumber\\
&=& h \left\{ \frac{1}{6}f(x_0) + \frac{2}{3}f(x_0+\frac{h}{2}) + \frac{1}{3}f(x_1) + \frac{2}{3}f(x_1+\frac{h}{2}) + \cdots \right. \nonumber\\
 && \left. + \frac{2}{3}f(x_{n-1}+\frac{h}{2}) + \frac{1}{6}f(x_n)    \right\}
\end{eqnarray}
这就是Simpson方法。
\begin{lstlisting}
/*
 * Simpson：辛普生积分
 * 给定积分区间[a,b]，区间分割为n份，被积函数为func，返回积分值
 */
double Simpson(double a, double b, int n, double (*func)(double)){

        double y1,y2,y,h;
        int i;
        h = (b-a)/n;//步长
        y1 = 0;
        for(i=0;i<n;i++){
                y1 += func(a+ (i+0.5)*h);
        }
        y2 = 0;
        for(i=1;i<n;i++){
                y2 += func(a+ i*h);
        }
        y = 4*y1 + 2*y2 + func(a) + func(b);
        y *= h/6;
        return y;
}
\end{lstlisting}

\subsection{误差分析}
利用泰勒展开，很容易推导出Newton-Cortes积分公式的误差。
为了方便叙述，我们记$m = x_i + h/2$。
我们写出$f(x)$在$m$附近的
\begin{eqnarray}
f(x) =&& f(m) + f'(m)(x-m) + \frac{f''(m)}{2} (x-m)^2  \nonumber\\ 
        &&+ \frac{ f^{(3)}(m)}{6}(x-m)^3 + \frac{ f^{(4)}(m)}{4!}(x-m)^4 + \cdots
\end{eqnarray}
其中 $f^{(n)}(m)$ 表示函数在 $x=m$ 处的 $n$ 阶导数。
所以有
\begin{eqnarray}
f(x_i) = && f(m) + f'(m)\frac{h}{2} + \frac{f''(m)}{2} \frac{h^2}{4}  \nonumber\\ 
        && + \frac{ f^{(3)}(m)}{6}\frac{h^3}{8} + \frac{ f^{(4)}(m)}{4!}\frac{h^4}{16} + \cdots \\
f(x_{i+1}) = && f(m) - f'(m)\frac{h}{2} + \frac{f''(m)}{2} \frac{h^2}{4}  \nonumber\\ 
        && - \frac{ f^{(3)}(m)}{6}\frac{h^3}{8} + \frac{ f^{(4)}(m)}{4!}\frac{h^4}{16} + \cdots \\
\end{eqnarray}

$f(x)$在区间$\left[x_i, x_{i+1}\right]$上的积分为
\begin{eqnarray}
\int^{x_{i+1}}_{x_i} f(x)dx &=& f(m)h  + \frac{f''(m)}{24} h^3 + \frac{ f^{(4)}(m)}{1920}h^5 + \cdots 
\label{equation:IntegralOn1Step}
\end{eqnarray}

将$f(x_i), f(m), f(x_{i+1})$的泰勒展开分别带入Newton-Cortes公式，得到它们的展开形式。
中值法
\begin{eqnarray}
\int^{x_{i+1}}_{x_i} f(x)dx &\approx& f(m)h,
\end{eqnarray}
梯形法
\begin{eqnarray}
\int^{x_{i+1}}_{x_i} f(x)dx &\approx& f(m)h + f''(m)h^3/8 + f^{(4)}(m)h^5/(4!2^4) + \cdots , \nonumber\\
\end{eqnarray}
辛普生方法
\begin{eqnarray}
\int^{x_{i+1}}_{x_i} f(x)dx &\approx& f(m)h + f''(m)h^3/24 + f^{(4)}(m)h^5/(1152) + \cdots , \nonumber\\
\end{eqnarray}
与公式（\ref{equation:IntegralOn1Step}）相比，可知中值法、梯形法的误差均为$O(h^3)$，辛普生方法的误差为$O(h^5)$。
但这是每一步长上的积分误差，所以总误差为：中值法、梯形法$O(h^2)$，辛普生方法$O(h^4)$。

\subsection{自适应算法}
当被积函数$f(x)$震荡的特别厉害的时候，我们需要比较小的步长$h$；如果被积函数$f(x)$非常平缓，较大的步长也可以得到不错的精度。
于是我们希望，程序根据函数的性质，自动选择合适的步长，进行积分。

我们以辛普生方法为例，上面推导出，它在$[x_i,x_{i+1}]$的误差
\begin{equation}
S^{[x_i,x_i+h]}_f - I^{[x_i,x_i+h]}_f = c h^5,
\end{equation}
其中$S_f$表示新普生方法得到的数值积分，c为常数(我们假定$f^{(4)}(x)$在$[x_i,x_{i+1}]$区间内变化很小。）

那么，如果我们将步长减半，则会得到
\begin{eqnarray}
&& S^{[x_i,x_i+h/2]}_f - I^{[x_i,x_i+h/2]}_f = c h^5/32, \nonumber\\
&& S^{[x_i+h/2,x_i+h]}_f - I^{[x_i+h/2,x_i+h]}_f = c h^5/32,
\end{eqnarray}
即
\begin{eqnarray}
S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - I^{[x_i,x_i+h]}_f = c h^5/16,
\end{eqnarray}
所以
\begin{eqnarray}
S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - S^{[x_i,x_i+h]}_f = - \frac{15}{16} c h^5,
\end{eqnarray}
即
\begin{eqnarray}
&& S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - S^{[x_i,x_i+h]}_f \nonumber\\
=&& - 15 (S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - I^{[x_i,x_i+h]}_f),
\end{eqnarray}
所以，如果我们要求$[x_i,x_{i+1}]$内的积分误差不超过$\epsilon$，则只需要求
\begin{eqnarray}
|S^{[x_i,x_i+h/2]}_f + S^{[x_i+h/2,x_i+h]}_f - S^{[x_i,x_i+h]}_f | < 15 \epsilon
\end{eqnarray}

我们可以写一个递归函数，实现自适应积分，代码如下
\begin{lstlisting}
/*
 * function Simpson applies the Simson's rule to calculate integration of func() over [a,b], breaking [a,b] into n small pieces.
 */
double Simpson(double a, double b, int n, double (*func)(double)){
        double h = (b-a)/n;

        double y = func(a) + func(b);
        for(int i=0;i<n;i++){
                y += 4*func(a+ (i+0.5)*h);
        }
        for(int i=1;i<n;i++){
                y += 2*func(a+ i*h);
        }
        y *= h/6;
        return y;
}
int count=0;
/*
 * SelfAdapt_Simpson uses self adapted algorithm to calculate integration of func() over [a,b].
 * It automatically breaks [a,b] into unequal pieces, according to the local situaion of func().
 * the error is controled to be less than epsilon
 */
double SelfAdapt_Simpson(double a, double b, double (*func)(double), double epsilon){

        count++;

        double y1 = Simpson(a, (a+b)/2, 1, func);
        double y2 = Simpson((a+b)/2, b, 1, func);
        double y3 = Simpson(a, b, 1, func);
        if( fabs(y1+y2-y3) < 15*epsilon ){
                return y3;
        }
        else{
                return SelfAdapt_Simpson(a, (a+b)/2, func, epsilon/2)
                        +
                        SelfAdapt_Simpson((a+b)/2, b, func, epsilon/2);
        }
}
\end{lstlisting}
其中，Simpson函数实现Simpson积分，而SelfAdapt\_Simpson函数实现自适应的Simpson积分。

我们可以做一个例子，比如做如下积分
\begin{equation}
\int^{0.1}_{10} \sin(1/x^2) dx,
\end{equation}
如图\ref{figure:integration_demon_self_adapt}所示，被积函数$\sin(1/x^2)$在图左侧振荡得非常厉害，在右侧却非常平缓。
所以在左侧需要很小的步长，在右侧稍大的步长即可。
\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{integration_demon_self_adapt.eps}
\caption{自适应辛普生积分：$\int^{10}_{0.1} \sin(1/x^2) dx$}
\label{figure:integration_demon_self_adapt}
\end{figure}
我们分别用两种方式做这个积分，一种是变步长自适应辛普生，另一种是同步长辛普生方法，得到的结果如表\ref{table:self-adapt}。
自适应积分方式，可以省下一个数量级的步数。
\begin{table}[h]
\centering
\begin{tabular}{cccc}
\hline\hline
\multicolumn{2}{c}{自适应} & \multicolumn{2}{c}{普通} \\
\hline 
步数 & 误差 & 步数 & 误差 \\ 
\hline 
1E4 & 4.71E-10 & 1E4 & -2.5E-6 \\ 
\hline 
 &  & 1E5 & -2.24E-10  \\
\hline 
\end{tabular} 
\caption{自适应辛普生积分：$\int^{10}_{0.1} \sin(1/x^2) dx$}
\label{table:self-adapt}
\end{table}

\section{高斯积分}

\subsection{多项式作为近似}

如果在$[a,b]$内用N阶多项式$P(x)$作为$f(x)$的近似，让$P(x)$经过$(x0=a,y0=f(a)), \cdots, (x_N=b, y_N=f(b))$一共N+1个点，则在$[a,b]$区间内任意一处，有
\begin{equation}
f(x) = P(x) + \frac{f^{(N+1)}(\epsilon)}{(N+1)!} (x-x_0)\cdots(x-x_N),
\end{equation}
积分以后误差为
\begin{eqnarray}
\epsilon_P &=& \int^b_a P(x)dx - \int^b_a f(x)dx        \nonumber\\
&=& \int^b_a \frac{f^{(N+1)}(\epsilon)}{(N+1)!} (x-x_0)\cdots(x-x_N) dx \nonumber\\
&\leq& \frac{f^{(N+1)}(\epsilon)}{(N+1)!} (b-a)^{N+2}   
\end{eqnarray}
所以区间$[x_i,x_{i+1}]$内积分的误差不超过$O(h^{N+2})$。
如果每一个小区间$[x_i,x_{i+1}]$内的积分都用N阶多项式作为近似，则整段积分的误差不超过$O(h^{N+1})$。
梯形法N=1，误差为$O(h^2)$，新普生方法N=2，所以误差为$O(h^4)$，小于$O(h^3)$。

所以，如果能用高阶的多项式作为函数的近似，一般来说，就会得到更小的误差。

\subsection{勒让德多项式}
为什么要讲这个？我们先卖一个关子，回顾一下勒让德多项式。
勒让德多项式是一类微分方程的解
\begin{equation}
C(1-x^2)P - m^2_l P + (1-x^2)\frac{d}{dx}((1-x^2)\frac{dP}{dx}) = 0.
\end{equation}
其中$m_l=0$时，解是勒让德多项式；$m_l\neq0$时，解是连带勒让德多项式，也有人叫做络和勒让德多项式，英文是Associate Legendre Polynomials。

勒让德多项式为
\begin{equation}
L_k (x) = \frac{1}{2^k k!} \frac{d^k}{dx^k} (x^2-1)^k, k=0,1,2,\cdots
\end{equation}
它具有正交性
\begin{equation}
\int^1_{-1} L_i(x) L_j(x) dx = \frac{2 \delta_{ij} }{2i+1},
\end{equation}
前面几项勒让德多项式为
\begin{eqnarray}
& L_0(x) = 1, \nonumber\\
& L_1(x) = x, \nonumber\\
& L_2(x) = \frac{1}{2}(3x^2-1), \nonumber\\
& \cdots \nonumber
\end{eqnarray}

\subsection{分解高阶多项式}
前面我们解释了，如果想要提高积分精度，可以用高阶的多项式作为函数的近似。
如果积分区域为$[-1,1]$，要求
\begin{equation}
\int^1_{-1} f(x) dx,
\end{equation}
我们想用$P_{2N-1}$，即2N-1阶多项式来作为$f(x)$的近似，我们可以将它分解为
\begin{eqnarray}
P_{2N-1}(x) = L_{N}(x) P_{N-1}(x) + Q_{N-1}(x),
\end{eqnarray}
这是可以做到的。
左侧有$2N$个参数，右侧$P_{N-1}(x),Q_{N-1}(x)$分别有$N$个参数，这个分解是唯一确定的。

我们注意到，$L_0(x), \cdots, L_{N-1}(x)$分别为$0,1,\cdots,N-1$阶多项式，所以$P_{N-1}(x)$一定可以表示为
\begin{equation}
P_{N-1}(x) = \sum^{N-1}_{i=0} a_i L_i(x),
\end{equation}
所以
\begin{equation}
\int^1_{-1} L_{N}(x) P_{N-1}(x) dx
= \sum^{N-1}_{i=0} a_i L_N (x) L_i(x) dx = 0,
\end{equation}

\subsection{高斯积分}
积分近似为
\begin{eqnarray}
\int^1_{-1} f(x)dx \approx \int^1_{-1} (L_{N}(x) P_{N-1}(x) + Q_{N-1}(x)) dx.
= \int^1_{-1} Q_{N-1}(x) dx
\end{eqnarray}

如果将$Q_{N-1}(x)$也用勒让德多项式表达，有
\begin{equation}
P_{N-1}(x) = \sum^{N-1}_{i=0} b_i L_i(x),
\end{equation}
而 $L_0(x)=1$，所以有
\begin{eqnarray}
\int^1_{-1} f(x)dx \approx \int^1_{-1} Q_{N-1}(x) dx　＝　\int^1_{-1} L_0(x) Q_{N-1}(x) dx
= 2 b_0,
\end{eqnarray}
所以只要我们知道了$b_0$的值，就得到了积分的近似值。

怎样求取$b_0$？我们能计算的是$f(x)$。。。
\begin{equation}
f(x) \approx L_{N}(x) P_{N-1}(x)　+ Q_{N-1}(x),
\end{equation}
$L_N(x)$在$(-1,1)$上有$N$个零点,我们把它们叫做$x_0, \cdots, x_{N-1}$，
则有
\begin{eqnarray}
f(x_i) = Q_{N-1}(x_i) = \sum^{N-1}_{j=0} b_j L_j(x_i),
\end{eqnarray}
$L_j(x_i), i,j=0,\cdots,N-1$构成一个方阵$L$，可以证明其行列式不为0，所以
\begin{equation}
b_0 = \sum^{N-1}_{i=0} f(x_i) L^{-1}_{i0},
\end{equation}
积分近似为
\begin{equation}
\int^1_{-1} f(x)dx \approx 2 \sum^{N-1}_{i=0} f(x_i) L^{-1}_{i0}
= \sum^{N-1}_{i=0} \omega_i f(x_i).
\end{equation}

$x_0,\cdots,x_{N-1}$是固定的，$2L^{-1}_{i0}$也是固定的，所以我们可以求得它们，对任何$f(x)$应用。
下表是$N-1=5$的例子
\begin{table} [h]
\centering
\begin{tabular}{ccc}
\hline \hline
i & $x_i$ & $\omega_i = 2 L^-1_{i0}$ \\ 
\hline 
0 & -0.90617 98459 38664 & 0.23692 68850 56189  \\  
1 & -0.53846 93101 05683  & 0.47862 86704 99366  \\  
2 & 0.00000 00000 00000 & 0.56888 88888 88889  \\  
3 & 0.53846 93101 05683  & 0.47862 86704 99366  \\  
4 & 0.90617 98459 38664 & 0.23692 68850 56189  \\ 
\hline \hline
\end{tabular} 
\caption{高斯勒让德积分：$L_5(x)$的零点，与相应的高斯积分权重。}
\label{tab:GaussLegendre5}
\end{table}

当然，积分区域不一定是$[-1,1]$，而是$[a,b]$，我们可以做如下变换
\begin{equation}
t = \frac{ x- (a+b)/2 }{ (b-a)/2 },
\end{equation}
则t的积分区域为$[-1,1]$,所以
\begin{eqnarray}
\int^b_{a} f(x)dx &=& \frac{b-a}{2} \int^1_{-1} f( \frac{a+b}{2} + \frac{b-a}{2}t ) dt 
\nonumber\\
&\approx& \frac{b-a}{2} \sum^{N-1}_{i=0} \omega_i f( \frac{a+b}{2} + \frac{b-a}{2}t_i ),
\end{eqnarray}
其中$t_i$为$L_N(x)$的零点，$\omega_i$为相应的权重，$N=5$时，如表\ref{tab:GaussLegendre5}所示。

\begin{lstlisting}
#include<iostream>
using namespace std;
#include<iomanip>

const double x_node[5]={
        -0.906179845938664,
        -0.538469310105683,
        0,
        0.538469310105683,
        0.906179845938664};
const double omega[5]={
        0.236926885056189,
        0.478628670499366,
        0.568888888888889,
        0.478628670499366,
        0.236926885056189};

double GaussLegendre5(double a, double b, double (*func)(double)){
        double x, y=0;
        for(int i=0;i<5;i++){
                x = x_node[i] * (b-a)/2 + (a+b)/2 ;
                y += omega[i] * func(x);
        }
        y *= (b-a)/2;
        return y;
}

double f(double x){
        return x*x*x*x*x;
}

int main(){

        double a=0;
        double b=1;
        cout<<"Integration of f(x) from "<<a<<" to "<<b<<" is "<<endl;
        for(int i=1;i<=12;i++){
                cout<<"n="<<i<<"\t ";
                n=i;
                cout<<setprecision(15)<<GaussLegendre5(a,b,f);
                cout<<"\t "<< GaussLegendre5(a,b,f) - 1.0/(i+1)<<endl;
        }
        return 0;
}
\end{lstlisting}


根据我们的理论分析，用高斯积分方法，相当于用$2N-1$阶多项式拟合被积函数，
\begin{eqnarray}
P_{2N-1}(x) = L_{N}(x) P_{N-1}(x) + Q_{N-1}(x),
\end{eqnarray}
在积分公式中，涉及到$L_N (x)$的$N$个零点，以及与之相关的$N$个权重。
\begin{eqnarray}
\int^b_{a} f(x)dx 
&\approx& \frac{b-a}{2} \sum^{N-1}_{i=0} \omega_i f( \frac{a+b}{2} + \frac{b-a}{2}t_i ),
\end{eqnarray}
上面我们写了$N=5$的代码，那么也就是说，我们用$9$阶多项式拟合被积函数。

那么我们得到一个结论：如果被积函数是多项式，阶数不大于9，我们的代码会非常精确；如果阶数大于9，理论上就有误差。
\begin{table}
\centering
\begin{tabular}{ccc}
\hline \hline
$f(x)$ & $\int^1_0 f(x)dx$数值解 & 误差 \\ 
\hline 
$x$     & 0.5  & -2.22044604925031e-16 \\ 
 
$x^2$ & 0.333333333333333 & -1.66533453693773e-16 \\ 
 
$x^3$ & 0.25 & -1.66533453693773e-16 \\ 
 
$x^4$ & 0.2 & -1.66533453693773e-16 \\ 
 
$x^5$ & 0.166666666666667 & -1.38777878078145e-16 \\ 
 
$x^6$ & 0.142857142857143 & -1.11022302462516e-16 \\ 
 
$x^7$ & 0.125 & -1.2490009027033e-16 \\ 
 
$x^8$ & 0.111111111111111 & -1.11022302462516e-16 \\ 
 
$x^9$ & 0.0999999999999999 & -1.11022302462516e-16 \\ 
 
$x^{10}$ & 0.0909076593600402 & -1.43154905070386e-06 \\ 

$x^{11}$ & 0.0833254598135549 & -7.87351977839246e-06 \\

$x^{12}$ & 0.0768984805856498 & -2.45963374271296e-05 \\
\hline\hline
\end{tabular}
\caption{5阶高斯积分：$\int^1_0 f(x) dx$，被积函数为多项式时的数值误差。}
\label{tab:Order-5-Gauss-Legendre}
\end{table}
表\ref{tab:Order-5-Gauss-Legendre}可以看到，计算结果与我们的理论预期是完全一致的，对于不超过9阶的多项式，5阶高斯积分的结果几乎等于理论值，只有一点舍入误差。

在实际应用中——我们再重复一遍——对于震荡比较厉害的被积函数，Newton-Cortes积分更加安全。
如果只取有限的几个点，进行高斯积分，是危险的。
这是因为，高斯积分的出发点是用高阶多项式拟合被积函数，而震荡很厉害的被积函数，可能需要非常高阶的多项式才能拟合好。。。

对于比较平缓的被积函数，用多项式拟合往往效果不错，所以可以使用高斯积分，只需取很少的几个点，可能非常高效。

当然，你也可以将积分区间分为很多份，然后在每个小区间里进行高斯积分。。。
